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PHOTOMETRIC REVERBERATION MAPPING OF THE BROAD EMISSION LINE REGION IN QUASARS 
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ABSTRACT 

A method is proposed for measuring the size of the broad emission line region (BLR) in quasars using 
broadband photometric data. A feasibility study, based on numerical simulations, points to the advantages and 
pitfalls associated with this approach. The method is applied to a subset of the Palomar-Green quasar sample 
for which independent BLR size measurements are available. An agreement is found between the results of the 
photometric method and the spectroscopic reverberation mapping technique. Implications for the measurement 
of BLR sizes and black hole masses for numerous quasars in the era of large surveys are discussed. 

Subject headings: galaxies: active — methods: data analysis — quasars: emission lines — surveys — tech- 
niques: photometric 



1. INTRODUCTION 

Black holes (BH) are physical entities that provide insight 
into the foundations of (quantum-) gravity theories. Such ob- 
jects are characterized by three basic properties: mass, angu- 
lar momentum, and charge. Among these properties the most 
fundamental is the mass, which characterizes all BHs. Angu- 
lar momentum is only relevant for rotating (Kerr) BHs while 
charge is less relevant for astrophysical BHs. 

Currently, the best evidence for the existence of BH is astro - 
physical. In particular, a broad range of BH masses is implied 
by astronomical observations: from solar mass BHs as the 
end products of massive stars to supermassive BHs (SMBH) 
at the centers of galaxi es with masses that may easily exc eed 
a million solar masses (lSalpeterlll96fl lLvnden-Belfl969l 

Perhaps the first evidence for the presence of SMBH at the 
centers of galaxies is associated with quasars and other types 
of active galactic nuclei (to which we shall henceforth refer to, 
generally, as quasars). In those objects, gas from the medium 
surrounding the SMBH, that is able to lose angular momen- 
tum and spiral inwards, eventually accretes onto it. As the 
gas accretes, its gravitational energy is effectively converted 
to radiation with most of the emitted flux originating from the 
immediate vicinity o f the SMBH, where the ga s reaches the 
highest temperatures (Shakura & Sunyaev 19731). 

The mass of SMBHs in small samples of low red- 
shift non-active galaxies has been me asured using, for ex- 
ample, maser- and stellar-kinematics ([Mivoshi et al.l 119951: 
iGenzelet al] 1 19971; iGebhardt et al.1120001) . For most low red- 
shift galaxies, SMBH mass estimation is, however, indirect 
and involves scaling relations that were derived for small sam- 
ples of obj ects for which dire ct mass measurements were pos- 
sible (Ma gorrian et al.lll998l) . 

As one attempts to estimate the mass of black holes at 
higher redshifts, various complications arise: it is not clear 
that locally determined scaling laws apply also to higher red- 
shift objects since galaxies and their environments evolve with 
cosmic time. Furthermore, high redshift objects discovered 
by flux-limited samples are often much brighter than their 
low-z counterparts, and it is yet to be established whether ex- 
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trapolated relations are at all relevant in that regime dNetzeri 
2003). For these reasons, direct SMBH mass determination 
of intermediate- and high-redshift objects is the most reliable 
way to further our understanding of SMBH-galaxy formation 
and co-evolution, and to reach a reliable census of black holes 
in the universe. 

Direct mass determination of SMBHs at high redshifts is 
currently limited to quasars, which are among the bright- 
est objects in the universe and are powered by the SMBH 
itself. While quasars certainly trace the galaxy popula- 
tion, they also present a particular epoch in a galaxy life- 
time, in which the black hole is accreting gas, and is at 
a stage of rapid growth. Therefore, better understand- 
ing of the SMBH demography in low- and high-z quasars 
has far-reaching impl ications for galaxy/ black-hole formation 
and cq - evolut i on dFerrarese et al] 120011 IShields et al 1 1 2003 ; 
l Osmerl |2004t iHopkins et all 12006b iKollmeier et all 2006 : 
Shields et all I2006t iGreene & Hoi l2006t IWoo et ail 12008 : 
iTrakhtenbrot & Netzerll2010l) . 

A reliable method to directly measure the mass of SMBHs 
in quasars involves the reverberation mapping of the broad 
emiss i on line region (BLR) in those objects (Bahc aTTet alJ 
[T97llPeterso"nl[l993l) . This region is k nown to be photo- 
ionize d by the central continuum source (Bald win & Netzer] 
119781) . and its emission, in the form of broad permitted emis- 
sion lines, is seen to react to continuum variations. Studies 
have shown that the BLR lies at distances, Rblr, that are 
much larger than the SMBH Schwarzschild radius, and also 
larger than the size of the continuum emitting region, i.e., the 
accretion disk. That being said, the BLR lies at small enough 
distances so that its kinematics is primarily set by the SMBH, 
which has a prominent contribution to the gravitational poten- 
tial on those scales. That this is the case may be deduced from 
the velocity dispersion of the emitted gas, as inferred from the 
width of the broad emission lines, which is typically of order 
a few x 10 3 km s _1 , as well as from the fact that, in luminous 
quasars, line variations lag behind relatively rapid continuum 
variations with typical time delays, fi ai > ~ /?bt,r Ic ~ 200 days, 
where c is the speed of light (Kaspi et al. 2000, and references 
therein). 

Technically, reverberation mapping involves multi-epoch 
spectroscopic observations so that continuum and emission 
line v ariations may be individually traced dDietrich et al.l 
[19931 iNetzer & Petersonl I1997L iKaspieTall 120001) . Contin- 
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uum and emission line light curves can then be used to in- 
fer the size of the br oad emission line region us ing a cross- 
correlation analysis (Bl andford & McKeei [19821 and §2.2). 
By measuring Rblr and the velocity dispersion of the gas (via 
spectral modeling of the emission line profile and by defining 
a measure of the velocity dispersion of the gas), the SMBH 
mass may be deduced, up to a geometrical factor of order 
unity (Peter son & Wa ndel 19991 lKrolikl2001). Mass determi- 
nations using this technique have been shown to be consistent 
with the results of other independent method s dGebhardt et al.l 
I2000t iFerrarese et"aT1l200U IPetersonll2010l). and are thought 
to be accurate to within 0.5 dex ([Krolik 200j]). 

Nowadays, reverberation mapping has a central role in de- 
termining SMBH masses in active galaxies and for shed- 
ding light on the physics of the central engine in these ob- 
jectfl Nevertheless, despite several decades of spectroscopic 
monitoring of quasars, reliable BLR sizes and correspond- 
ing SMBH mass measurements ha ve been determined for 
only a few dozens o f low- z objects (iGaskell & Sparkdl 19861: 
Korat kar & Gaskelll [19911; iDietrich et al.l Il993h llCaspi et al] 
20001: IPeterson et al.ll2004t IKaspi et alj 120051; IPeterson et alj 
20051 iBentz et all 120091: iPetersonl l2010t iBarth et al.1 1201 ll) . 
Using the results of reverberation mapping, various scal- 
ing relations were found and quantified including, for ex- 
ample, the R m n (or SMBH mass) vs . quasar luminos- 
ity relation (iKoratkar & Gaskelll 119911: IKaspi et . al.l 120001: 
IPeterson et al Il2004t IKaspi et alJl2007l;lBentz et alj|2009l) . the 
SMBH mass vs. the qua sar's optical spectral properties re- 
lation dGreene & Holl200"5l) . a nd the black h ole-host galaxy 
bulge mass relation in quasars (lWandellll999h . 

Despite being derived from small samples of nearby 
objects, the aforementioned relations are often extrapo- 
lated to reflect on the quasar population as a whole 
dVester gaard 2002; Vestergaard & Osmer 2 0091). sometimes 
with little justifica t ion dNetzerl 120031: IKaspi et alj 120071: 
iMcGill et alJ 120081: IZhang et alJ l2008h . Potential biases 
arising from uncertain extrapolations to high redshifts, or 
to different sub-classes of quasars, may crucially affect 
our understanding of black hole growth and structure for- 
mation over cosmic time d Kauffmann & Hae hneltl 120001; 
Netzerl 120031: ICollin & Kawaguchil 12004 iPelupessv et all 
20071: iTrakhten brot & Netzerll2010l) . Clearly, an efficient way 
to implement the reverberation mapping technique is highly 
desirable as it would alleviate many of the uncertainties cur- 
rently plaguing the field. 

Motivated by upcoming photometric surveys that will cover 
a broad range of wavelengths and will regularly monitor a fair 
fraction of the sky with good photometric accuracy (e.g., the 
Large Synoptic Survey Telescope, LSST), we aim to provide 
proof-of-concept that photometric reverberation mapping of 
the BLR in quasars is feasible. In a nut shell, instead of 
spectrally separating line and continuum light curves from 
multi-epoch spectroscopic observations, we take advantage of 
the different variability properties of continuum and line pro- 
cesses and separate them at the light curve level. In this work 
we have no intention to address the full scope of issues associ- 
ated with time-series analysis but rather hope to trigger further 

3 It is worth noting that, in addition to SMBH mass measurement, the 
structure and geometry of the BLR may be probed by mean s of velocity- 
resolved reverberation mapping of the broad emission lines {Korista et all 
[T993 ; IWanders et al.llT995 ; Kollatschny 2003; B entz et alj|2010f) . Such u> 
vestigations cannot be carried out using broadband photometric reverberation 
techniques, which do not resolve the emission lines, and will not be discussed 
here any further. 
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Figure 1. A realization of the broadband photometric light-curves of a typ- 
ical quasar (see text). Flux variations in the X (black) and Y (red) bands 
show an overall similarity due to the large contribution of highly correlated 
continuum emission processes to both bands. The response of an emission 
line to continuum variations in the X-band is shown in dashed blue lines 
(tf = 200 days was assumed). The Y-band light-curve includes a 13% con- 
tribution to the total flux from a broad emission line (solid blue line). A case 
where the noise level is negligible (5/ (/) = 10~ 3 ) was assumed for clarity. 

observational and theoretical work in broadband reverberation 
mapping of quasar BLRs. 

This paper is organized as follows: section 2 presents the 
photometric approach to reverberation mapping using ana- 
lytic means and by resorting to simulations of quasar light 
curves. The reliability of the photometric method for line- 
to-continuum time lag measurement is investigated by means 
of simulations in §3. In §4 we apply the method to a subset 
of the Palomar-Green (PG) quasar sample (Schm idt & Greenl 
1983), fo r which previous measurements of the BLR size are 
available (Kaspi et al. 2000). Summary follows in §5. 

2. THE FORMALISM 

Here we present the algorithm for reverberation mapping 
the broad emission line regions in quasars using broadband 
photometric means. We first describe our model for quasar 
light curves (§2.1) and then analyze them using new statistical 
methods that are developed in §2.2. 

2.1. A model for the light-curve of quasars 

To demonstrate the feasibility of photometric reverberation 
mapping, we first resort to simulations. Specifically, we rely 
on current knowledge of the typical spectrum of quasars, and 
the variability properties of such objects, to construct a model 
for the light curves in different spectral bands. 

2.1.1. Continuum emission 

Quasar light-curves are reminiscent of red noise, following 
a power spectru m P(co) ~ (O a ((0 is the angular frequency) 
with a ~ —2 dGiveon et al.l fl999, and references therein). 
Here we model the pure continuum light-curve in some band 
X as a sum of discrete Fourier components, each with a ran- 
dom phase 0,-: 

/i(0 = E A ( w ') sin ( fi) ' f +^)- (d 

;=] 

cc /2 

The amplitude of each Fourier mode, A(fijf) °= » ( and t is 
tima3 We define the normalized variability measure, % — 

4 We note that, realistically, the power at each frequency is distributed 
around the above-defined powerlaw, forming an effective envelope with some 
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Figure 2. The emission line transfer functions used in this work to model the 
broad lines' response to continuum variations (see text). While a gaussian 
transfer function (thick line) is relatively well localized in the ti me domain, 
the transfer function which is reminiscent of those discussed in M aoz et aT] 
(1991) extends to much longer times. The latter transfer function naturally re- 
sults in long-range correlations between line and continuum processes (whose 
effects are discussed in §3.2.5), and roughly characterizes the reaction of a 
geometrically thick BLR that is isotropically encompassing the central con- 
tinuum source. We note that both transfer functions are normalized such that 
the line-to-continuum time-delay is ~ ?j' ae . 

oj - S 2 / (/} [c.f. iKaspi et all (120001) who express this in 

per-cent], where (/} is the mean flux of the light-curve, ay is 
its standard deviation, and 8 is the mean measurement uncer- 
tainty of the lig ht curve. For the Palom ar-Green (PG) sam- 
ple of quasars (Sch midt & Greenlll983l) . t he continuum nor- 
malized variability measure % c ~ 0.1 -0.2 (IKaspi etalJ2000h 
and, typically, 8 > 0.01 (see table 1). A light-curve realiza- 
tion, with 5 ^ ay and % c ~ 0.2, is shown in Fig. Q] 

Quasar emission occurs over a broad range of photon en- 
ergies, from radio to hard X-ray wavelengths and originates 
from spatially distinct regions in the quasar. In the opti- 
cal bands, there is s ome evidence for r ed continua to lag 
behind blue continua (lKroliketai]ll991t ICollier et aTlll998h 
ISereeev et al.l2005l:lBachevl2009l) . Without loss of generality, 
we associate the T-band with the redder parts of the spectrum 
so that its instantaneous continuum flux may be expressed as 
a function of the flux in the X-band, 



j*{t)=f x *xf 



dxf x (x)y c (t-T), 



(2) 



where \j/ c (8t) is the continuum transfer function^. Presently, 
\j/ c is poorly constrained by observations and only the 
inter-band continuum time delay, f lag , is loosely determined 
(Bachev 2009, and refer ences therein). Nevertheless, so long 
as Rblr/c ^ f lag , the particular form of \j/ c is immaterial to 
our discussion (see §3.2.3). For simplicity, we shall therefore 
take \j/ c to be a gaussian with a mean of ? lag and a standard 

deviation f( a( ,/2. 

2.1.2. Line emission 

Emission by photo-ionized (BLR) gas is naturally sensitive 
to the continuum level of the ionizing source. In particular, 

characteristic width. We ignore this additional complication here, which in- 
troduces yet another degree of freedom in the model but does not qualitatively 
change the results. Real light curves are analyzed in §4. 

5 All transfer functions, y/, satisfy the normalization condition 
f+:d{8t) V (8t) = 1. 



for typical BLR densities and distances from the SMBH, the 
response time of the photoionized gas to continuum variations 
is much shorter than the light travel time across the region, 
and so may be considered instantaneous. Therefore, the flux 
in some line i, which is emitted by the BLR, can be written, 
in the linear approximation, as 



f!=fx*v4, 



(3) 



where the transfer function y/- depends, essentially, on the 
geometry of the BLR region a s seen from the direction of 
the observer dHorne et al.H2.Q04l and references therein). Our 
present knowledge of y/'(5f) for the quasar popu l ation as a 
whole, is, however, rather limited [see Goa d et al.l (11993b for 
a few relevant forms for x// 1 ]. In what follows we shall con- 
sider two forms for the transfer function: a) a gaussian kernel 
with a mean t! and a standard deviation fjL/4, and b) a trans- 
fer function with a flat core in the range < 8t < 1.36?/ , 
with a following linear decline to zero at St = 2.72f ]ag (such 
a transfer function also results in a time-lag being ~ t[ ). 
The latter transfer function is moti vated by the ob servations 
of the Seyfert 1 galaxy NGC 415 1 dMaoz et all 19911) and cor- 
responds to a case where a geometrically thick, spherical shell 
of broad emission line gas isotropically encompasses the con- 
tinuum emitting region. The different transfer functions de- 
fined above will henceforth be termed (a) the gaussian and (b) 
the thick-shell BLR models. The transfer functions are de- 
picted in figure [2] The main qualitative difference between 
the two transfer functions, as far as light curve behavior is 
concerned, is that while the gaussian one is relatively local- 
ized in time, the thick shell BLR model is much broader, and 
results in longer range correlations between the light curves 
of the two bands. Unless otherwise stated, our results are pre- 
sented for the gaussian model. Cases for which we find a 
qualitative difference between the two transfer functions are 
explicitly treated (e.g., §3.2.5). 

A simulated light-curve of an emission line that is driven by 
continuum variations (assuming f lag = 200 days) is shown in 
figure Q] (blue curves). Time-delay and smearing effects, with 
respect to the continuum light-curve, are clearly visible. Here 
too, we define a normalized variability measur e for the line, 
Xi, and consider Xi — 0.75 Xc (Kasp i et alJ2000l see their table 
5). The sum of continuum and lines' emission to the total flux 
in band Y is therefore given by 



(4) 



where f Y is the total flux in the emission lines. While a sim- 
ilar expression applies, in principle, also for the total flux in 
bandX, unless otherwise stated, we shall assume that fx = f x , 
i.e., that emission line contribution to band X is either negli- 
gible or that the line variability associated with it operates on 
timescales that are much too short to be resolved by the ca- 
dence of the experiment, or too long to be detectable in the 
time-series. For simplicity, we shall limit our discussion to a 
single emission line contributing to the T-band (this restric- 
tion is relaxed in §3.2.4). 

The contribution of the emission line(s) to the total flux in 
the F-band requires knowledge of the equivalent width of the 
broad emission l ine, W [typically < lOOA in the optical band 
for strong lines (IVanden Berk et all 1200 lb l. the transmission 
curve of the broadband filter, as well as the redshift of the 
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Figure 3. Statistical estimators for evaluating the line to continuum time delay. Left: The ACFs for the X and 7 bands are shown (same color coding applies 
as that in Fig. [TJ, as well as the CCF of the two bands (black curve). Note the excess power in the ACF of the Y-band, as well as in the CCF of the two bands, 
relative to the ACF of the Z-band (gray hatched region). This results from relatively long-range correlations between line and continuum variations, which occur 
on BLR light travel times. The excess power is the sought after signal (see text). Right: The statistical estimators, ^ca, ^AA, which quantify the excess power 
due to line emission, are consistent with f/ c (see legend). The position of the peaks at positive time-lags (the only ones that preserve causality) is identified 
with the line-to-continuum time-delay. Evidently, all methods give consistent predictions for the time-lag, which is also in agreement with the input value, t t 
(dash-dotted black line). 

quasar. For example, for filters with sharp features in their 
response function, even small redshift changes could lead 
to substantial differences in the relative contribution of the 
emission line to the total flux in the band, rj. For a fiat fil- 
ter response function over a wavelength range Wj, (typically, 
W h ~ 10 3 A for broadband filters), 77 = W /Wt- (A more ac- 
curate treatment of the line contribution to the broadband flux 
for the case of Johnson-Cousins filters is given in §4.) Exam- 
ple light-curves for the X and Y bands, assuming 77 = 0.133, 
are shown in figure [TJ In the following calculations, we shall 
work with normalized light-curves that have a zero mean and 
a unit variance, i.e., / — > (/ — (/))/ Of. 

2.2. Measuring the Line to Continuum Time Delay 

The spectroscopic approach to measuring the size of the 
BLR cross-correlates pure line and continuum light-curves so 
that the position of the peak (or the center of mass) of the 
cross-correlation function (CCFffl E,i r (8t) = fy * f\, marks 
the line to continuum time-delay (Peterson 1988, and refer- 
ences therein). However, line and continuum light curves can- 
not be disentangled using pure broadband photometric means 
and the data include only their combined signal. In fact, light- 
curves in different bands are rather similar due to the dom- 
inant contribution of highly correlated continuum emission 
processes (Fig. [TJ. Therefore, a naive analysis in which the 
CCF between the light curves of different bands is calculated 
would reflect more on f-fL than on fi (see Fig. 0). 

Seeking to measure the time-delay between line and con- 
tinuum processes using broadband photometric means, it is 
important to realize that instead of spectrally distinguishing 
between line and continuum light curves, it is possible, un- 
der certain circumstances, to separate these processes at the 
light curve level. In particular, the line light curve, f Y = 
fy — f§ — fy — fx since f Y ~ This is justified for quasars 



6 Evaluation of t he cros s - and auto- correlation functions is done here 
via the definition of Welsh ( 1999, see his equation 6) for the local cross- 
correlation fjrnctioji^_Unless otherw ise stated, we use the interpolated CCF 
method of Gaskell & Sparge fi986l) . 

7 This is expected to be true at least for adjacent bands since the signal is 
largely due to continuum emission processes, such as black-body radiation, 
which cover a broad wavelength range. Interestingly, this broad similarity 



since the continuum auto-correlation timescale [> 10 2 days 
dGiveon et al.l 19991) . possibly reflecting on viscous timescales 
in the accretion flow] is much longer than f[ ag [being of order 
days (ISergeev et al.l l2005) and resulting from irradiation pro- 
cesses in the accretion disk, which take place on light travel 
time scales]. Furthermore, as fx = ./£, f Y — fy ~ fx, and we 
obtain, 

~ &a{X,Y) ee {f Y -fx) *fx = fr*fx-fx *fx, (5) 

where the rightmost term is just the difference between the 
CCF of the light-curves and the auto-correlation function 
(ACF) of the light curve with the negligible line contribu- 
tion (see §3.2 for cases in which emission lines contribute 
also to fx). Graphically, this difference is simply the shaded 
region shown in the left panel of figure [3] which arises due 
to correlated line and continuum processes adding power on 
?j -scales. Using these statistics, and given a rather non- 
restrictive set of assumptions (see more in §3), it seems, 
in principle, possible to recover the line-to-continuum time- 
delay using pure broadband photometric data. In particular, 
for the example shown in figure [3] ^ca peaks at the correct 
time-lag, and is in agreement with the time-lag estimate as 
obtained using the spectroscopic method, i.e., using £,i c . 

The estimator defined by equation [5] provides a good ap- 
proximation for t,i c in cases where the emission line contribu- 
tion to the band is sub-dominant (see below). For the more 
general case (e.g., when narrow band filters are concerned), 
it is possible to define a somewhat more generalized version 
where 4ca = fy * fx — (1 — * fx, which is applicable 
for cases in which the emission line contribution to the Y- 
band is considerable or even dominanfl Nevertheless, with 
this definition, 77 needs to be observationally known, which is 
not always the case, especially when large (photometric) red- 
shift uncertainties are concerned, the filter response function 

of light curves in different bands seems to be true also for the (hidden) far- 
UV and optical bands given the success of the spectroscopic reverberation 
mapping technique in measuring line to continuum time-delays. The effects 
of inter-band time-delay and time-smearing of the continuum signal may ac- 
count for the observe d structure function differences between adjacent bands 
I MacLeod et al. 2010), and are further discussed in §3.2.3. 

8 In this case, lim^-ji ££ A = fy c . For r\ -C 1, as appropriate for broadband 
photometry, & A ~ £ C A- 
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has sharp features (IChelouche et al.l201 ll) . or the quasar spec- 
trum, if available, is noisy. Focusing on broadband photomet- 
ric light curves in this work, the contribution of emission lines 
to the band is, typically, of order a few per cent (see §4). We 
therefore restrict the following discussion to £,qa as defined 
by equation|5] 

It is possible to define an additional statistical estimator to 
measure line to continuum time-delays using broadband pho- 
tometry. Specifically, by using the fact that, for broadband 
filters and to zeroth order approximation, fx — fy (see Fig. 
[TJ, we may approximate ^ca by 

$AA(X,Y)=fY*fr-fx*fx- (6) 

This estimator is just the difference between the ACFs of the 
light-curves in each band. A more intuitive way to under- 
stand the above definition may be gained by looking at fig- 
ure [3] where the ACF of filter Y shows more power on ? lag 
timescales due to relatively long-term correlations between 
line and continuum processes. In contrast, the ACF of filter 
X, being devoid of line emission, does not show a correspond- 
ing power-excess. We note that both definitions of £, are quasi 
anti-symmetric with respect to X and Y exchange. 
We find that, as far as ? lag measurements are concerned, both 

estimators, £,qa and <^aa> give. m principle, consistent results 
(see Fig. |3). While £,qa has the advantage of being a more 
reliable estimator of §/ c (see more in §3), <^aa has the advan- 
tage of not containing a cross-correlation term of the form 
fy * fx- This property allows one to evaluate ^aa even in 
cases where non-contemporaneous light curves are available 
so long as a) both light curves are adequately sampling the 
line transfer function, and b) quasar variability results from 
a stationary process (see §3.2.5). Moreover, if the variability 
properties for an ensemble of objects (of some given prop- 
erty) are statistically similar, and the auto-correlation terms, 
fy * fy and fx *fx, are well determined in the statistical sense 
(see §3.1.1), then it is not required for the same objects to be 
observed in both bands. This opens up the possibility of mon- 
itoring different parts of the sky in each band, and using the 
acquired data to statistically measure the mean line to con- 
tinuum time-delay. That being said, large ensembles may be 
required to achieve a reliable measurement of f' ag since ^aa 
is, in fact, an approximation to an approximation of §/ c , and 
its results are expected to be somewhat less reliable than both 
<^CA and §/ c (this is further discussed in §3). 

A note is in order concerning the maximal value of £,^(8t) 
and |cA(5f). Unlike the spectroscopic technique where 
%lc(8t) is, in fact, the Pearson correlation coefficient and its 
value measures how well continuum and line processes are 
correlated at any 8t, for the photometric estimators defined 
here, the peak value depends also on the relative contribution 
of the emission line(s) to the broad band flux. The latter de- 
pends not only on the properties of the quasar, but also on 
the transmission curve of the broadband filter, the variability 
measure of the line, and the redshift of the object. Therefore, 
the significance of the result is not directly implied by the 
value of the peak in <!;aa and E,qa, and quantifying it requires 
a different approach, as we shall discuss in §§3,4. 

3. SIMULATIONS 

While the reliability of the spectroscopic reverberation 
mappi ng technique has been a ssessed in numerous works 
re.g.. iGaskell & Petersonl (1 1987b and references therein], we 
have yet to demonstrate it for the photometric method. In this 



section we take a purely theoretical approach and use a suite 
of Monte Carlo simulations to test the method and quantify its 
strengths and weaknesses. 

3.1. Sampling & Measurement Noise 

All light curves are finite with respect to their duration time, 
f to t, as well as the sampling interval, f sam . This fact alone intro- 
duces potentially substantial "noise" in the cross-correlation 
analyses (whether photometric or spectroscopic). The effect 
of finite sampling is seen in Fig. |4] (left panels) where the 
light curves of individual objects are analyzed, showing that 
<^ca, <^aa formally peak at values that could be significantly 
removed from f lag . This implies that even in the ideal case, 
where measurement noise is negligible, the deduced time-lag 
might be significantly offset from the true value. While sim- 
ilar situations occur also with spectroscopic data, our simula- 
tions indicate that the photometric estimators are, statistically, 
more prone to such effects. 

In addition to finite sampling there is measurement noise, 
which reflects on the ability of any statistical estimator to re- 
cover the true time-lag with good accuracy. On an individual 
case basis, the effect of noise is seen to give rise to more er- 
ratic <^ca and ^aa behavior, potentially leading to substantial 
deviations in the position of the peak from the true time-lag. 
In fact, for the particular case shown (Fig. |4), £,qa has a kink 
at ? lag and peaks only at ~ 3f lag . In this particular case, i^aa 
happens to be less affected by noise, and formally peaks at 
~ 1 • 15f lag due to a narrow feature being the result of noise. 

Sampling and measurement noise are manifested as small 
scale fluctuations in all statistical estimators. In particular, in 
some cases (see Figs. H] & [5]), large amplitude fluctuations 
are observed, which peak at long timescales, <>?/f lag > 1, and 
their value could exceed the value of the true peak (i.e., that 
which is associated with, and peaks at, ~ t lag ) . Clearly, an 
algorithm that would aid to properly identify the relevant sig- 
nal, would be of advantage. To this end we use a physically- 
motivated algorithm related to the geometry of the BLR: the 
fact that the BLR encompasses the continuum emitting region 
means that the width of its line transfer function should be of 
order the time-lag. Therefore, one expects the relevant sig- 
nal in i^cA) ^aa to peak at f lag with a width, A(8t) ~ ? lag . In 
particular, narrower features with A(5f)/f la <C 1 are likely to 
be merely the result of noise or sampling, and although they 
could reach large amplitudes, they are of little relevance to the 
time-lag measurement. 

The aforementioned properties of the line transfer func- 
tion motivate us to consider the following algorithm to aid in 
the identification of the time-lag: we resample £,qa and <^aa 
onto a uniform logarithmic grid, which preserves resolution 
in A(8t)/8t units. We then apply a running mean algorithm 
such that the boxcar size corresponds to a fixed range in log- 
arithmic scale. This procedure smooths out narrow features 
having A(8t)/8t <C 1. Examples for smoothed statistical es- 
timators are shown for two levels of noise in figure Clearly, 
narrow features at large 8t are significantly suppressed with- 
out diminishing the power of similar A(8t ) features centered 
at smaller 8t. It should be noted, however, that using this 
technique does not assure that the correct time-lag is recov- 
ered (see above). Furthermore, in some cases, multiple peaks 
may be found and additional methods are required to assess 
their significance (see §4). 



6 Chelouche & Daniel 

1 10 100 1000 




-2 2 -2 2 ( f 2 2 -2 2 

Figure 4. Averages of ijcA (blue) ,^aa (red) as well as for (black) for ensembles of statistically equivalent quasars (the number of objects in each ensemble 
is indicated above each column). Upper (lower) panels correspond to a noise-level 8/ {/) = (0.03). In all simulations the total duration of the light curve is 
4f| ag , and visits are f/ ag /50 apart. The relative contribution of the emission line to the filter flux, rj ~ 0.13. Clearly, finite sampling of the light curves, even in 
the absence of measurement noise, could result in the peaks being offset with respect to ( la This is true for individual objects as well as for small samples. 
Generally, however, upon averaging over many objects, the mean estimators peak at St ~ t( (right panels). When measurement noise is present, larger ensemble 
are required for the peak of the {£ca) , (£aa) and (§; c ) to coincide with / lag to good accuracy (e.g., compare the top and bottom panels in the second column 
from the left). Overall, the photometric and the spectroscopic estimators give consistent results, which is especially true for large ensembles. 
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Figure 5. Photometric time-lag measurements for individual quasars having 
noisy data (5/ (/) = 0.03/0.01 in the top/bottom panel). Calculated ^ca's 
(blue curves) show multiple peaks, some of which are narrow and peak at 
St > tfag, potentially with large amplitudes (note the feature at St/t[ ig ~ 3). 
Smoothing of the signal using a running mean algorithm over a uniform 
logarithmic grid suppresses narrow features and highlights the peak associ- 
ated with the time lag (cyan curves; see text). Red arrows mark the peak at 
'lag measured' which, in this example, deviates only slightly from the input lag 
(dashed vertical line at 8/t^ = 1). 

Testing for the accuracy of time-lag measurements using 
noisy data, we have conducted uniformly sampled light curves 
simulations of a quasar (100 realizations were used having 
800 points in each light curve, with a total light curve dura- 
tion, ?tot, satisfying ?totA' ag = 5 and 77 ~ 0.13) at several levels 
of measurement noise. For each realization, ^ca was calcu- 
lated and St where it peaks was identified with f| measured and 
logged. The resulting time-lag statistics are shown in figure 
[6] The distribution is clearly non-gaussian and shows a peak 
at the input time lag with a tail extending to long time-lags. 
Not unexpectedly, the dispersion in time-lag measurements 
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Figure 6. Measured time-lag distributions for a sample of 100 quasars (hav- 
ing 800 points in each light curve and with ftot/'/ ag = 5; 7] ~ 0.13 was as- 
sumed), and for several noise levels (see legend). Clearly, for all cases, the 
distributions peak at the input time lag (dashed line), demonstrating the fea- 
sibility of the photometric approach for time-lag measurement. As expected, 
noisier light curves result in broader time lag distributions. The distribu- 
tions show a tail extending to large / lag me asured Al'ag~ vanies an ^ we clearly 
non-gaussian. 

increases with the noise level. For the particular examples 
shown, taking 8/ (f) = 0.01 [a noise level typica l of pho- 
tometric light curves of quasars (Kaspi 'eTaIll2000l) l. ~ 60% 
of the measurements will result in measured being within 

±10% of the true lag. For noisier data with 8/ (/) =0.1, only 
~ 30% of all realizations result in a time lag within ±10% of 
the true time lag. These calculations demonstrate that photo- 
metric reverberation mapping in the presence of measurement 
noise typical of photometrically monitored quasars is indeed 
feasible, and that the measurement is, in principle, reliable. 
When noisier data are concerned, the uncertainty on the time- 
lag is larger. 
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Figure 7. Relative measurement uncertainties on t{ measured- Left: One standard deviation of the measured time-lag distribution functions (i.e., the uncertainty 
on 'lag measured) as obtained by measuring the time delay for individual realizations in a large ensemble of objects (having the same BLR size and sampling 
characteristics as those of figure[5]. One hundred simulations were used to calculate each data point (6 equally log-spaced points were calculated in each curve 
and a second degree polynomial fit to the results is shown). Noise level 8/ (/) = 0.01 (0.1) is shown in solid (dashed) lines, and the usual color coding scheme 
applies (^ca/^aa/^/c results in blue/red/black curves). Clearly, better signal-to-noise is obtained by either reducing the measurement noise or by better sampling 
the light curve. Generally, £ca is the more reliable photometric estimator, while both photometric methods are less accurate than the spectroscopic one (£/<■). 
Right: The uncertainty on the mean time-lag deduced for an ensemble of objects using ijcA, for two noise levels (5/ (/) = 0.01. 0.1 are shown in solid and 
dashed curves, respectively, and smoothing was applied). The total observing effort (number of points for the entire sample) was held fixed. The time-lag was 
deduced in two ways: (1) by measuring the time-lag for individual objects with the £,ch estimator and taking the mean result (thick lines), and by measuring 
the position of the peak of (£ca) (thin lines). The uncertainty on the mean time-lag reported for the latter quantity was deduced from multiple simulations of 
statistically equivalent ensembles. The results show that, for the same observing effort, monitoring more objects with lower cadence (so long as the line transfer 
function is properly sampled) results in more accurate time-lag measurements. Also, measuring the peak in the mean statistical estimator provides more accurate 
results than measuring the time-lag for individual objects [o"((^ca)) < CT ('iag measured)!- 



To better quantify the ability of the photometric method to 
uncover the time-lag under various observing strategies, we 
repeated the above analysis for two noise levels (8/ {/) = 
0.01,0.1) while varying the number of points in the light 
curves (i.e., observing visits). More specifically, we main- 
tained a fixed ftotAi' ag = 5 and varied the sampling period, 
fsam (we made sure that f lag /4am > 5 so that the BLR transfer 
function is never under-sampled). Monte Carlo simulations 
were carried out to obtain reliable statistics, and we report 
the one standard deviations of the measured time-lag distribu- 
tions, <7(fi' a g,measured)' in figure0(also shown, for comparison, 
are the corresponding ^-statistics). As expected, better sam- 
pled light curves result in a lower dispersion in fj ag measured for 

all statistical methods. Also, unsurprisingly, <^ca is the more 
reliable photometric estimator for the time-lag, and both pho- 
tometric estimators are worse than §; c . In particular, quasar 
light curves having ^100 points in each filter with photomet- 
ric errors of order 1% result in a ^cA-deduced time-lag dis- 
persion of order ~ 35 %. Similar quality light curves can b e 
found in the literature dGiveon et al.lll999tlKaspi et alJ uOOO), 
and are analyzed in §4. 

Consider the following question: having finite resources, 
is it better to observe fewer objects with better cadence or 
more objects having sparser light curves? [We assume that 
the cadence is adequate for sampling the line transfer func- 
tion (see more in §3.2), and that the same signal-to-noise is 
reached in both cases.] The answer to this question is given 
in figure [7] (right panel) where the total number of visits was 
kept fixed yet the number of objects in the ensemble varied 
(correspondingly, the number of visits per object changed to 
maintain their product constant). Interestingly, for the same 
observing effort, it is advantageous to observe more objects 
with worse cadence. Using this strategy, one is less sensi- 
tive to unfavorable variability patterns of particular quasars. 
For the case considered here, the time-lag uncertainty drops 



by a factor ~ 3 by increasing the sample size from one to 
ten objects (the uncertainty for the case of ensembles was de- 
termined from Monte Carlo simulations of many statistically 
similar ensembles). Similar results are obtained also for <^aa 
(not shown). 

3.1.1. Ensemble averages 

In the upcoming era of large surveys, statistical information 
may be gathered for many objects. For example, the LSST 
is expected to provide photometric light curves for > 10 6 
quasars having a broad range of luminosities and redshifts. 
With such large sets of data at one's disposal, sub-samples of 
quasars may be considered having a prescribed set of prop- 
erties (e.g., luminosity and redshift), and the line to contin- 
uum time-lag measurement may be quantified for the ensem- 
bles. In this case, it is possible to measure the characteris- 
tic line-to-continuum time-delay in two ways: by measuring 
the peak of 4ca(5?) or §aa(5?) on a case by case basis and 
then constructing time-lag distributions to measure their mo- 
ments, or by defining the ensemble averages of the statistical 
estimators. Specifically, we define (£,cA(8t)) , (§aa(5?)) as 
the arithmetic averages^ of the corresponding statistical esti- 
mators at every 8t. (Alternative definitions for the mean do 
not seem to be particularly advantageous and are not consid- 
ered here.) The average statistical estimators for ensembles of 
varying sizes are shown in figure |4] As expected, averaging 
leads to a more robust measurement of the peak, as the various 
<jjcA(5f)'s "constructively-interfere" at ~ f' ag . The ensemble 
size required to reach a time-lag measurement of a given ac- 
curacy naturally depends on the noise level. 

We find that measuring the time-lag from the peak of (<^ca) 
(or (^aa)) is more accurate than measuring the time-lags for 

9 We similarly define, for comparison purposes, (£,i c (8t)). These defini- 
tions were, in fact, previously used to compute the curves shown in figure 
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Figure 8. The effect of scatter in the broad line region properties on (arbi- 
trarily scaled) (£rj A (5f)). A large ensemble of objects (typically 100-1000) 
was used, whose time-lag distribution follows a modified Gaussian (see text), 
which is characterized by the mean time-lag, fi ag and its standard deviation, 
O'C'lag)- f totAiag — 4 was assumed. For c(fi ag ) /?i ag < 1, the distribution is 
nearly Gaussian, {^cA(St)) is, likewise, symmetric, and its peak is at fi ag . 
For CT(Zlag)/f]ag ^ 1. the time-lag distribution is, by definition, modified (see 
text), and" the effective mean time lag is > ?] ag (black squares for each curve). 
This, and the fact that, for some objects, ~ f tot (hence the line transfer 
function is poorly sampled) result in an asymmetric (£ca(5?)), which only 
roughly peaks at the effective mean time lag. 

individual objects in the ensemble, and then quantifying the 
moments of the distribution (Fig. FFP\ This is particularly 
important when noisy data are concerned, in which case reli- 
able time-lag measurements for individual objects are difficult 
if not impossible to obtain. To see this, we note the rising un- 
certainty on the deduced time-lag, a(t[ mea sured)> as tne num ~ 
ber of visits per object is reduced (and despite the fact that a 
correspondingly larger sample of objects is considered; see 
the right panel of Fig . |7]for a case in which 8/ (/) =0.1). In 
contrast, when working with (^ca) -statistics, the uncertainty 
monotonically declines so long as the sampling of the transfer 
function is adequate. 

The above analysis highlights the importance of using large 
samples of objects (even with less than ideal sampling) in ad- 
dition to studying individual objects having relatively well- 
sampled light curve. In particular, by using large enough sam- 
ples of quasars, the effects of noise, and even sparse sampling, 
may be overcome. 

Current studies suggest that there is a scatter o f 30%-40% 
in the BLR size vs. qu asar luminosity relation (Ben tz et al.l 
20091 lKaspiet1d]|2000h . To assess the effect of such a scat- 
ter on (£ca) (similar results apply also for (§aa})> we aver- 
aged the results for an ensemble of simulated objects whose 
time-lags follow a Gaussian distribution with a mean t\ ag and 
standard deviation a(t\ a g). We modify this Gaussian such that 
the input lag is always positive definite: for <7(fi a2 )Aiag <C 1 
the distribution is Gaussian while for <7(/i a g)/fi ag ~ 1 there 
is an excess of systems with zero lag. (^cA(St)}, is shown 
in figure H] evidently, the average statistical estimator peaks 
at the mean of the distribution so long as C7(fi ag )/fi ag < 0.8. 
For larger ratios, the distribution considerably deviates from 
Gaussian and its mean is > ?i ag . In this case, (|cA(5f)) still 
peaks around the mean but its shape is highly asymmetric. 

10 The uncertainty on the time-lag, as deduced from (^ca), was calcu- 
lated by simulating numerous statistically-equivalent ensembles of objects, 
calculating {^ca)> and measuring the time-lag from its peak. A time-lag 
distribution is obtained with the uncertainty given by its standard deviation, 



This is partly due to the time-lag distribution, but also so due 
the fact that the transfer functions with the largest fj ag 's are not 

adequately sampled by the light curve (?t t/? lag = 4 was as- 
sumed). To conclude, photometric reverberation should yield 
accurate mean BLR sizes for a sample of quasars, given the 
observed scatter in the BLR properties. 

3.2. Systematics 

We identified several potentially important biases in the 
method, which we shall now quantify using a suite of Monte 
Carlo simulations. We work in the limit of large ensembles 
so that the effect of sampling and measurement noise are neg- 
ligible (typically, the results of several hundreds of realiza- 
tions are averaged over and the line to continuum time delay 
is identified with the peak of (§ca(^)) and (§aa(^))). 

3.2.1. Measurement noise 

Biases are present in the time-lag measurement from noisy 
data. These are negligible for noise levels 8/ (f) <r\ but may 
surface otherwise. Specifically, for a (Gaussian) noise level 
of ~ 50%, and a relative line contribution to the band of 10%, 
both photometric estimators over-estimate f| by as much as 

20% with ^aa being more biased than ^ca (Fig- HJ- This bias 
could therefore be present in situations where large and noisy 
ensembles are considered, or when weak emission lines are 
targeted. Quantifying the bias requires good understanding of 
the noise properties. 

A related bias exists when the noise level in the two bands 
differs substantially. More specifically, when \8x/ (fx) — 
Sy/ ify) I ~ i? {8x/8y is the noise level in the X-/Y- bands), 
our simulations show that fi agjme asured ma y differ from ?i ag by 
as much as 25%. Whether an over-estimate or an under- 
estimate of fi ag is obtained depends on the particular noise 
properties, as mentioned above. Observationally, one should 
therefore aim to work with data having comparable noise lev- 
els for the different bands to obtain a bias-free result. As noted 
above, for relatively strong lines or low noise levels, such bi- 
ases are negligible. 

Varying seeing conditions leading to aperture losses, as well 
as improperly accounting for stellar variability in the field 
of the quasar, could lead to correlated errors between the 
bands. This in turn could lead to artificial zero-lag pea ks in 
the correlation function of the light curves (<Welshlll999|) . and 
seve ral solutions have been devised to overcome this prob- 
lem (lAlexandeiill997t[Zu et al.ll201 ll and references therein). 
Nevertheless, the formalism developed here subtracts off, by 
construction, the contribution of similar processes (e.g., cor- 
related noise) to both light curves. In fact, as far as the algo- 
rithm is concerned, the continuum contribution to the 7-band 
is merely an effective correlated noise to be corrected for, 
which it does. Therefore, photometric reverberation mapping, 
as proposed here, is considerably less susceptible to the ef- 
fects of correlated noise than the spectroscopic method. That 
being said, if the noise between the bands correlates over suf- 
ficiently long (e.g., ^>day) timescales, a secure detection of 
the line-to-continuum time delay may be more challenging 
(for an effectively similar case see §3.2.3). 

3.2.2. Sampling 

Sampling is rarely uniform. In fact, seasonal gaps could 
be of the order of the time lag in some objects (Kaspi et al. 
2000). The full treatment of sampling related issues is be- 
yond the scope of this paper, and we restrict our discussion to 
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Figure 9. Biases in measuring line to continuum time delays as calculated for 
the spectroscopic (abscissa) and photometric methods (ordinate). Dashed line 
is for a 1 : 1 ratio. Upper-left panel shows the effect of noise with the largest 
data point [blue (red) reflecting on ^ca-(^aa-) deduced measurements] cor- 
responding to 8/ (f) = 0.512 with subsequently smaller size symbols corre- 
sponding to values smaller by a factor 2 (the same symbol-size value applies 
to all panels). Upper-right panel shows the bias as a function of < Bap and im- 
plies consistency between the spectroscopic and photometric results (largest 
symbol is for %apA]' ag = 2). Lower-left panel quantifies the bias for the X- 
band leading (circles) or trailing (squares) the F-band. Largest symbol is for 
I'faglAliio — 0'3. Lower-right panel shows the bias in the presence of a sub- 
ordinate emission line whose time-lag is t[ /2, which is contributing to the 
Y (circles) and to the X (squares) filter. The largest symbol corresponds to 
J] s /j] — 0.8 (see text). 

light curve gaps (real data are treated in §4). To check for po- 
tential systematics, we resampled our simulated light curves 
with equal "on" and "off observing periods with durations 
fgap- We identify a bias when f gap ~ f' . Nevertheless, both 
the spectroscopic and the photometric methods are equally af- 
fected by it so that their results remain consistent. As our 
aim is to focus on cases in which the photometric approach 
leads to new or different biases, we do not further quantify 
this bia{3 

3.2.3. Inter-band continuum time delays 

Thus far, tf <C t[ was assumed. While the spectroscopic 
method is unaffected by inter-band continuum delays (the 
continuum in the immediate vicinity of the emission line is 
considered), our definition of the continuum is that which cor- 
responds to the filter having the least contribution of emission 
lines to its flux, i.e., the X-band. 

For a large set of simulation runs (obtained by varying both 
rj and ?i ag Ai ag ; not shown), we find that biases become non- 
negligible for rj -1 |?f ag |A' ag > 1 (ff ag > implies that contin- 
uum emission in the F-band lags behind the X-band). Here we 
consider particular cases in which | ? iag I / ^lag — ^-3 (i) — 0.13 
is assumed), and find that fi ag |Ai ag > 0.1 results in substan- 
tial, > 20% bias (Fig. |9]l for the following reason: although 
modest values of l^glAiag are considered, the contribution of 
continuum processes to both light curves by far exceeds that 
of the emission line. Specifically, when the continuum in the 
F-band lags behind the X-band, the deduced time lag reflects 
more on f lag than on t' . When the X-band lags behind the 

F-band, and f^A/ag is significant, the peak in both statisti- 

1 ' See §4 for the usage of discrete correlation functions, which are less 
susceptible to sampling issues. 



cal estimators is strongly suppressed and time-lag measure- 
ments become less reliable. Realistically, inter-band contin- 
uum time-delays seem to be of order 10~ 2 f lag dBachevll2009l) . 
and it is therefore likely that biases of the type discussed here 
would be negligible for most strong broad emission lines. 

3.2.4. Multiple emission lines 

Realistic quasar spectra could results in more than one 
emission line contributing to the F-band, or even different 
emission lines contributing to both bands. While such compli- 
cations are irrelevant for the spectroscopic method (note the 
small statistical scatter around unity in Fig. |9), the photomet- 
ric method is only sensitive to the combined contribution of 
all emission processes in some waveband. Here we study the 
biases introduced by including the contribution of a weaker 
secondary emission line to either of the bands whose time lag, 
t' ^t 1 

r lag,s ^'lag- 

We find potentially considerable biases when the relative 
contribution of the weaker line to the flux, rj s satisfies rj s / rj < 
1. However, the magnitude of the bias is also affected by 
('lag s — f iag) / f iag- We ^° not attem pt t0 fully cover the param- 
eter space in this work, and focus on particular cases. Figure|9] 
shows an example with t[ s /f lag = 0.5 where we find > 20% 
biases for rj s /rj > 0.5: when the weaker line contributes to 
the F-band, f lag measured < t{ while the opposite occurs when 
the weaker line contributes to the X-band. Comparing to a 
case in which f lags = ? / ag / 1 and r\ s /r\ = 0.8, we find sim- 
ilar trends but an overall smaller ( < 5%) bias. Assuming 
instead f lag s = 2f ]ag and r\ s /r\ = 0.8, we find a bias of order 
20% (10%) when the weaker line contributes to the Y- (X-) 
band. While the magnitude of the bias is clearly not solely 
determined by rj /rj s , we find that the bias is, generally, small 
for ri s /ri <C 1. Clearly, if f/ ag s < f sam or s > t to t then the 
time-series is insensitive to the presence of the secondary line, 
and the bias would be negligible regardless of rj s /rj. 

Further quantifying the bias as a function of the line proper- 
ties (e.g., transfer functions and time-delays), filter response 
functions, quasar redshifts, and observing patterns is beyond 
the scope of this paper, and should be carried out on a case by 
case basis. 

3.2.5. Total light curve duration 

It is well known that the ratio of the total duration of the 
light curve to the time-lag, ?totAi' ag ' needs to be greater than 
some value to be able to adequately sample the transfer func- 
tion and reliably measure the time lag using spectroscopic 
means. Our calculations demonstrate that the same holds also 
for the photometric method. In particular, for the gaussian 
(thick shell) transfer function, reliable measurements require 
thaUtotALg £ 3 (9); see figureE)] 

Differences between the photometric and the spectroscopic 
approach arise for ftotA' ag ^ 1 ; while the spectroscopic 
method gradually converges to q (assuming sufficient data 
and adequate sampling; see Fig. [101 . the photometrically de- 
duced time lag, ? lag measure( j, develops a bias with respect to 
f lag with increasing ftotAi ag - Generally, the bias is comparable 
(but not identical) for (^ca) an d (^aa). an d over-estimate the 
true lag. We quantify the bias for the two line transfer func- 
tions, y/ 1 , defined above (§2.1.2), and find it to be different 
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Figure 10. Left: the bias in the time-lag measurement as a function of the observable quantity h ot /t[ measure( j (^ca/^aa/^/c results are shown in blue/red/black 
curves). For all statistical estimators there is a minimum 'tot Ai' ae measured below which the sampling is inadequate. Above this threshold, the spectroscopic method 



converges to the input time lag (e.g., for a gaussian transfer function, t[ r 



^lagforftotAlag,: 



method wherein tL 



lag .measured / lag 



/t( is a rising function of 'totAiagr 



ig, measured 



> 3). In contrast, there is a bias using either photometric 



j, and is more pronounced for geometrically thick BLR configurations. Right: a particular 



: 2fj (dark-black solid line). 



example for (ijcA(<5')), calculated for an ensemble of 100 statistically identical objects with 'tot Ai' ag = 50, resulting in a peak at St '. 
When the light curves of both bands are de-trended by a polynomial of degree d, the bias is gradually removed. For large enough rf-values, substantial power 
on relevant timescales is removed, thereby affecting the significance of the results (see legend). When the quasar power-spectrum is truncated at frequencies 
CO < (Omi,, [here (27r/(D mm )Atot — 0.13], the light curve appears to be stationary, and the bias disappears (dash-dotted dark blue line). 

in each case: broader Xj/ 1 that extend over larger 8t intervals, 



lead to more biased results. For example, while the thick shell 
BLR model shows a factor ~ 2 bias for ?totAi ag .measured = 50, 
the corresponding bias for the Gaussian BLR is only ~ 40% 
(see Fig. flOl . The bias is seen to grow logarithmically with 
f totA/ag measured ' anc ' so would t> e most prominent when ana- 
lyzing long time-series of low-luminosity quasars. Consider- 
ing, for example, ~ 10 years of LSST data for a low-z quasar 
with a thick shell BLR whose size is ~ 100 light-days. In this 
case, ftot/f^ ~ 30, which would result in ^measured A/ ag ~ 2. 
Not correcting for this bias would lead to an over-estimation 
of the black hole mass, by a similar factor, with considerable 
implications for black hole formation and evolution. 

Consider the bias in (^ca) : as its value is sensitive to the 
line transfer function, its origin may be traced to the fy * fx 
term i n Eq. [5] which is akin to the bias discussed in Welsh 
( 1999). In a nutshell, with the quasar's power spectrum being 
red, most of the power is associated with the lowest frequency 
modes. Nevertheless, those are the same modes that are not 
adequately sampled by any finite time series. As such, the 
data appear to originate in a non-stationary process where, for 
example, the mean flux changes substantially between the ex- 
treme ends of the light curve. This and the fact that emission 
lines' contribution to the flux is always additive results in a 
positive bias. More extended line transfer functions, resulting 
in longer term correlations between the light curves, are there- 
fore more prone to such biases. With (^aa(^O) — (4cA(5f))> 
similar biases are also encountered there, that trace back to 
the fy * fy term. 

To verify the source of the bias, we repeated the calculations 
using a truncated form for the power-spectrum where all fre- 
quencies below a minimum frequency, (D m i n , have no power 
associated with them (we have experimented with 27r/fi ag 3> 
Omin 3> 27r/ftot)- The resulting (^ca(5?)) is indeed unbiased 
and peaks at fj ag (not shown). For intermediate cases where, 
for example, the po wer spectrum flattens at low frequencies 
dCzernv et al.lll999l) . the results will be less biased but may 
not be entirely bias-free. 



Having identified the source of the bias, it is now pos- 
sible to correct for it using one of several methods which 
we term: de-trending, sub-sampling, and bootstrapping. De- 
trending has been shown t o be useful t o correct for biases in 
the spectroscopic method (Welsh 1999). Basically, one filters 
out long-term trends in the light curve by fitting (done here 
in the least-squares' sense) and subtracting off a polynomial 
of some degree, d, from the data. As a polynomial of de- 
gree d has d zeros, subtracting a higher degree polynomial 
from the data results in a suppressed variance at frequencies 
CO < d(2n/ttot)- Conversely, the highest degree polynomial 
that can be subtracted from the data, whilst not suppressing 
the reverberation signal, is d < [ f tot/fiagJ- The effect of de- 
trending is shown in figure [To] for a range of li-values, and 
assuming ftotAiaa = 50. While subtracting a d = 2 polynomial 
already decreases the bias substantially, a bias-free result is re- 
covered by de-trending with a d = 20 (= 0.4f tot /f la J polyno- 
mial. Further increasing d suppresses the sought after signal 
to insignificant levels by removing modes with frequencies 
~ 27r/f lag . The dependence of the remaining bias on the de- 
trending polynomial degree, d, is shown in the inset of figure 

023 

An alternative method to remove the bias is to analyze con- 
tinuous sub-samples of the time series, thereby reducing the 
effective t tot and decreasing the bias. This method can be im- 
plemented in cases where one is not data-limited. 

The third method for removing the bias involves bootstrap- 
ping: given an observed continuum light curve (i.e., fx), and 
an assumed iff 1 , the light curve of the line-rich band is sim- 
ulated. Upon sampling it using the cadence of the real data, 
an artificial fy light curve is obtained. By calculating the sta- 
tistical estimators for the semi-artificial data, it is possible to 
compare f lag measured to the input lag, and correct for the bias. 
Conversely, should the bias be known by other means, infer- 
ences concerning the nature of the line transfer function may 
be obtained. 

4. DATA ANALYSIS 
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Table 1 

Photometric properties of the PG sample of quasars'"' 



Object 




AL A (5100A) 




B 






R 




net 


net 


Ha lag 1 *' 


Predicted lag 


Photo, lag 


Name 




(10 44 ergs _1 ) 


pts 


Xb 


(8b) 


pts 


Xr 


(Sr) 


wo/Fe 


w/Fe 


(days) 


(days) 


(days) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(ii) 


(12) 


(13) 


(14) 


PG 0026+129 


0.142 


7.0±1.0 


70 


0.19 


0.02 


71 


0.14 


0.02 


0.03 


0.07 


132+g 


147 ± 25 


- 


PG 0052+251 


0.155 


6.5 ± 1.1 


75 


0.26 


0.03 


75 


0.20 


0.02 


0.04 


0.08 


2H±|| 


141 ±27 


- 


PG 0804+761 


0.100 


6.6±1.2 


83 


0.17 


0.01 


88 


0.14 


0.01 


0.05 


0.07 


193+20 


136 ±27 


200 ±100 


PG 0844+349 


0.064 


1.72 ±0.17 


65 


0.10 


0.02 


66 


0.10 


0.02 


0.08 


0.06 


Jv 16 


51±6 


500 ±200 


PG 0953+414 


0.239 


11.9±1.6 


60 


0.14 


0.03 


60 


0.12 


0.01 


0.06 


0.11 


1 87+^ W 


231 ± 39 


- 


PG 1211+143 


0.085 


4.93 ±0.80 


25 


0.16 


0.01 


24 


0.13 


0.01 


0.06 


0.07 


U6tl 


109 ±20 


- 


PG 1226+414 


0.158 


64.4 ±7.7 


45 


0.12 


0.01 


45 


0.10 


0.01 


0.04 


0.08 


5 14+65 


703 ±135 


- 


PG 1229+204 


0.064 


0.94 ±0.10 


44 


0.16 


0.03 


45 


0.16 


0.03 


0.08 


0.06 


-,+39 


34±4 


N/A 


PG 1307+085 


0.155 


5.27 ±0.52 


30 


0.14 


0.01 


30 


0.10 


0.01 


0.04 


0.08 


1/y -145 


122± 16 




PG 1351+640 


0.087 


4.38 ±0.43 


35 


0.11 


0.01 


35 


0.10 


0.01 


0.06 


0.07 


247+jf 


101 ± 13 




PG 141 1+442 


0.089 


3.25 ±0.28 


29 


0.10 


0.01 


29 


0.07 


0.01 


0.06 


0.07 


103+1° 


82±9 




PG 1426+015 


0.086 


4.09 ±0.63 


26 


0.16 


0.01 


25 


0.13 


0.02 


0.06 


0.07 


90+« 


96 ±17 




PG 1613+658 


0.129 


6.96 ±0.87 


66 


0.11 


0.01 


64 


0.10 


0.02 


0.04 


0.07 


«3 


144 ±22 


300 ±130 


PG 1617+175 


0.114 


2.37 ±0.41 


56 


0.19 


0.02 


56 


0.14 


0.02 


0.04 


0.07 


1 1 1+31 
1 -37 


67 ±12 




PG 1700+518 


0.292 


27.1 ±1.9 


53 


0.07 


0.02 


53 


0.07 


0.01 


0.05 


0.11 


11 4+246 (c) 
235 


428 ±61 


300 ±160 


PG 1704+608 


0.371 


35.6 ±5.2 


38 


0.13 


0.01 


40 


0.11 


0.01 


0.02 


0.05 


4^7+252 (c) 


550 ±109 


400 ±200 


PG 2130+099 


0.061 


2.16 ±0.2 


78 


0.10 


0.02 


79 


0.08 


0.02 


0.08 


0.06 


237+g 


60±7 


240 ±100 



Columns: (1) Object name, (2) redshift, (3) monochromatic luminosity in units of 10 44 erg s , (4)-(6) properties of the photometric light curve of the B-band: 
number of points, the normalized variability measure, and the mean measurement uncertainty. (7)-(9) having the same meaning as (4)-(6) but for the R band. (10) 
the net contribution of emission lines to the bands excluding the iron blends. (1 1) like (10) but including the iron emission complexes. (12) the spectroscopically 
m easured time lag fo r the Ha line from Kaspi et al. 1 2000). (13) The obsen>ed lag for the Balmer emission lines as predicted by the BLR size-luminosity relation 
of Kaspi et al. 1 2000, see their Eq. 6). (14) The time lag measured in this work using broadband light curves. 

(a) Highlighted table rows indicate objects for which a statistically significant signal is detected (see §4.3). 

(b) Values correspond to centroid observed-frame values (Kaspi et al. 2000, see their table 6). 

(c) Data for Ha are not available. Time lag corresponds to the H/3 line. 



Having gained some understanding of the advantages and 
limitations of broadband photometric reverberation mapping, 
we now wish to test the method using real data. Unlike the 
models described in §3, real data are often characterized by 
uneven sampling, gaps, and a varying signal to noise across 
the light curve. Here we consider a subset of the Palomar- 
Green (PG) sample of objectsEI (ISchmidt & Greenll 19831) . for 
which Johnson-Cousins B and R photometry is publicly avail- 
ably dMaoz et all 1 1991 iGiveon et alJ[l99l . and indepen- 
dent BLR size measurements are available throu gh spectro- 
scopic reverberation mapping (Kas pi et"aTl 12000). The prop- 
erties of the sample are given in table 1. 

While the PG sample provides a natural testbed for our pur- 
pose, it is far from being ideal: the typical number of points 
per light curve is small, ~ 50, and the photometric errors non- 
negligible, < 2% (c.f. Fig. [7]). Other photometric data sets 
are, in principle, available (Geha et al. 2003; Meusinger et al. 
1201 ll) but do not yet have BLR size measurements w i th whi ch 
our findings may be compared [see lChelouche et all (1201 1|) 1. 

4.1. Preliminaries 

We first identify line-rich and line-poor bands given 
the typical quasar spectrum ( Vanden Berk et al.l 1200 11) . the 
spectroscopically- determined redshift for the PG sample 
( Kaspi et al. 2000), and the relevant filter throughput curves 
(Fig. ITTb . In particular, we need to identify the filter with the 
larger relative contribution of emission lines to the variance 

12 We choose to not focus our attention on active galactic nuclei (i.e., low 
luminosity quasars) since, in this case, great care is needed to make sure that 
varying seeing conditions do not introduce additional noise whose character- 
istics are poorly understood. In this case, specific data reduction techniques, 
such as aperture photometry, may be of advantage. 

13 http://wise-obs.tau.ac.il/PG.html 



of the light curve. This is, however, difficult to estimate, since 
the variability properties of only a few lines, in only a few 
objects, are sufficiently well understood. Instead, we adopt a 
simplified approach and assume that the contribution of emis- 
sion lines to the variance is proportional to their relative con- 
tribution to the flux. This is a fair approximation if, for exam- 
ple, the relative flux variations in all e mission lines are com- 
parable. This is supported by theory (Kaspi & Netzerj [i999l) 
and by the fact that the fractional variability of th e Balmer 
lines is independent of their flux ( Kaspi et al. 2000, see their 
table 5). 

To properly account for the contribution of emission lines 
to the flux in some band, spectral decomposition of all promi- 
nent line, line blends, and continuum emission processes is 
required. To this end, we considered t he high signal-to- 
noise (S/N) composite quasar spectrum of Vande nBerk et al.l 
(120011) . and identified all prominent emission lines and blends 
in the relevant wavebands. To fit for continuum emis- 
sion, we identified line-free regions ove r the entire spec- 
tral range from rest-UV to the optical (Kaspi etaTI 120001; 
Vanden Be rk et al.l2001l) . and fitted a double powerlaw model 
(see Fig. [TTJ, which accounts for the changing slope at optical 
wavelengths. The stronger Balmer emission lines were fitted 
with a two Gaussian model, which accounts for their rela- 
tively narrow line cores and extended wings. Weaker lines 
were fitted with a single Gaussian model. Iron blend emis- 
sion templates, as derived fo r IZw I (IVeron-Cettv et al.ll2004t 
IVestereaard & W ilkes 2001), were convolved with a Gaus- 
sian kernel and fitted to the data. Independently adjusting the 
width and normalization for each Gaussian kernel (individual 
lines and blends), a reasonable agreement was sought. We 
model the B almer continuum using the analytic expression of 
iGrandil dl982l) . We emphasize that our model is by no means 
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Figure 11. Emission line and continuum contribution to the Johnson-Cousins broadband filters. Left: the composite quasar spectrum (Vanden Berk et al. 2001, 
assumed to be at z = 0.1; see table [TJ is shown together with the B- and i?-filter response functions. Spectral decomposition into prominent emission lines 
(magenta) and iron line blends (gray curve) is shown. A double powerlaw continuum emission model is also included (dashed and dash-dotted magenta lines), 
and the "best-fit" model is shown in green. Right: the net emission line contribution to the broadband light curves as a function of the quasar redshift (dashed 
black lines; gray hatched patterns indicate the uncertainties assuming a 10% uncertainty on the flux of individual lines). The net contribution of a few particular 
lines is shown in colored curves. The sum of all lines excluding iron (light grey) and including iron (dark grey) are shown, indicating that in either case the R-filter 
is relatively line-rich for z < 0. 1 objects with the Balmer emission line contribution being the domi nant one. Typical redshift intervals ove r which line-rich filter 
identifications change are greater than the typical uncertainty associated with photometric redshifts (Weinstein et al. 2004; Wu & Jia 2010). 

(RSS). In a nutshell, the first algorithm adds noise accord- 
ing to the measurement errors, and repeats the analysis many 
times to estimate the uncertainty associated with the results, 
while the latter method resamples the light curves (while los- 
ing ~ 37% of the data), to reach the same goal. Quite often 
both tests are combined to form the FR-RSS method, which 
is our method of choice in the present analysis. When con- 
sidering th e DCF approach, th e z-transformed DCF algorithm 
(ZDCF) by Alexander] (11997b often produces the best results, 
and is the approach adopted here. While all the aforemen- 
tioned algorithms have been thoroughly tested with respect to 
spectroscopic reverberation mapping, their adequacy for the 
photometric reverberation mapping method is yet to be exam- 
ined. 

Our implementatio n of the ICCF (FR and RSS algorithms 
alike) is standard dPeterson et alj Il998). Nevertheless, be- 
cause photometric reverberation mapping requires the sub- 
traction of two cross-correlation functions, and to avoid in- 
terpolations of any kincQ our implementation of the ZDCF 
algorithm requires that both the fy * fx and fx * fx terms in 
equation [5] are evaluated over an identical time grid, and that 
the bin center is identified with the time tag for that bin, i.e. 
8t. Similar considerations apply also when evaluating <^aa- 
Furthermore, we require that the number of independent pairs 
contributing to ea ch bin in the correlation functions is > 1 1 
(Alexander 1997). For simplicity, we use equally spaced bins 
to evaluate ^ca an d <^aa> an d repeat the an alysis many time s 
by randomly selecting independent pairs (Alexander 1997). 
This provides us with a distribution of <^ca an d <^aa values at 
each time bin, from which the mean and its uncertainty (asso- 
ciated with one standard deviation) are evaluated. 

Spectroscopic reverberation mapping often employs two 
distinct approaches to estimate the time lag: by measuring the 
position of the peak of the correlation function, or by mea- 
suring its centroid. While the results are comparable, they 
are not identical and many works have discussed the pros and 
cons associated with each approach dKaspi et al J 120001 and 



physical, nor uniquely-determined, and that its sole purpose is 
to allow for the qualitative estimate of the emission line con- 
tribution to the flux in each band. The decomposed spectrum 
is shown in figure [TT] 

The line-rich band, as a function of the quasar redshift, 
for the Johnson-Cousins filter scheme, is shown in figure [TT] 
Clearly, depending on the object's redshift, either the B or the 
R filter may be identified with the line-rich band. Typically, 
a redshift accuracy of 0.2 is required to securely identify the 
line-rich band, which is within reach of photo metric redshift 
determination techniques (Richards et al. 2001). We find that, 
for z < 0.5 quasars, the R band may be identified with the line- 
rich band (i.e., with the T-band, using the nomenclature of 
§2.1.2). This result is rather insensitive to whether or not the 
optical iron blends contribute to the variance on the relevant 
timescales. For z < 0.25 objects, the flux and, by assumption, 
the variance, are largely dominated by the Balmer lines. Un- 
certainties of order 10% in the flux of individual lines and line 
blends do not significantly change this result. 

4.2. Methods 

Aiming to analyze individual objects in the PG sample hav- 
ing a marginally adequate number of points in their light 
curves (thereby resulting in a large uncertainty on the deter- 
mined lag; Fig. |7), we outline the formalism used below to 
assess the significance of our results. 

4.2.1. Cross-correlation analysis 

In analyzing real data for individual objects in the PG sam- 
ple (table 1), it is important to verify that the detected sig- 
nal is real and is not some artifact of non-even sampling or 
the particular light curve observed. There are several meth- 
ods in the literature for cross-correlating realistic data sets: 
the in terpolated cross-correlation function (Gaskell & Sparke 
119861 ICCF; used in S3) an d the discrete correlation func- 
tion fedelson & Krolik|[T98l DCF). While the ICCF may be 
somewhat more sensitive to low-level sign als, it is also more 
prone to sampling artifacts (lPetersonll993l) . To verify the sig- 
nificance of the ICCF signal, two algorithms are often used to 
estimate the uncertainty in the deduced correlation function: 
the flux randomization (FR) and the random subset selection 



14 Although we do not use this method here, we note that due to the corre- 
lated nature of the correlation function, interpolation at the correlation func- 
tion level is likely to provide results which are less sensitive to sampling 
artifacts than the ICCF method. 
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Figure 12. Light curve analysis for PG 0804+76 1 . Left: The B and R light curves show significant non-monotonic variability over a few hundreds days timescales 
(see table 1). Right: the calculated statistical estimators, £ca (blue shades) and i^aa (red shades). Also shown is the renormalized (black points) based on 
the spectroscopic data of Kaspi et al. (2000). Two analysis methods are used to calculate the photometric estimators: ICCF (solid lines and shades) and ZDCF 
(discrete points), which lead to consistent results. The uncertainty on each estimator was calculated using 100 Monte Carlo simulations: the FR-RSS method was 
used to estimate the error on the ICCF, and a random independent pair selection was used to estimate the error on the ZDCF (see text). All estimators peak at 
St ~ 200 days. The result obtained using <^ca has a higher statistical significance, and both photometric estimators are less significant than the spectroscopic one. 

references therein). Similar algorithms are implemented here. M ' 

Estimating the uncertainty on the measured time-lag is done 
in the following way: we use the set of Monte Carlo simula- 
tions described above and calculate £,qa, <^aa f° r eacn real- 
ization. The time-lag is identified either with the peak or with 
the centroid of the statistical estimators^ 5 !, and a distribution 
of time-lags is obtained. The uncertainty on the time-lag is 
associated with one standard deviation of that distribution. 

4.2.2. Controlling systematics 

In addition to the above algorithms, aimed at testing the ro- 
bustness of the signal, we provide an additional test to identify 
systematic effects of the type discussed in §3.2.5 and to assist 
in screening against problematic data sets that may lead to 
spurious signals. To this end we use the bootstrapping method 
discussed in §3.2.5, verifying that the mean contribution of 
the simulated emission line to the flux of the line-rich band is 
that given in Fig. [TT] and that the variance o f the line light 
curve is consistent with the Kaspi et al. (2000, see their table 
5) results. 

4.3. Results 

Here we outline the results for individual objects enlisted in 
table 1: we first discuss quasars for which a significant signal 
is detected 1 ^, and then cases where the analysis does not result 
in a discernible peak in either estimator. Our treatment of the 
ensemble as a whole is given in §4.3.9. 

4.3.1. PG0804+ 761 

The B and R light curves for PG 0804+761 (henceforth 
termed PG0804), at z = 0.1, are shown in figure [T21 where 
large non-monotonic flux variations are e vident. This objec t 
has the best sampled light curves in the iKaspi etafl (2000) 
sample, with a mean uncertainty on the flux measurement 

15 Here we use the smoothing algorithm defined in §3.1, to avoid identify- 
ing spurious peaks with the physically meaningful peak. 

16 The significance criterion used here requires a > Iff signal (with o" 
being the one standard deviation uncertainty on the measured value) detected 
in both statistical estimators, and for which the ICCF and ZDCF analysis 
yield consistent results. 



tributions for PG0804, obtained using sets of 100 Monte Carlo simulations 
(standard color scheme is used). ICCF (ZDCF) results are shown in the upper 
(lower) panels. When using the ZDCF approach, two different bin sizes were 
assumed (solid and dashed lines; red/blue shades correspond to £aa/^caX 
and consistent results are obtained in both cases. The mean lag and its uncer- 
tainty (assumed to be one standard deviation) corresponding to each distri- 
bution are denoted in each panel. Centroid statistics results in slightly larger 
time-lags due to the non-symmetric shape of the statistical estimators around 
the peak (Fig. f5J. While the time-lags inferred by all statistical estimators and 
methods are statistically consistent, the uncertainties inferred by the ZDCF 
approach are relatively small, and possibly less reliable (see text). 

in both bands being of order 1 % (with maximum uncertain- 
ties on individual measurements being ~ 2%). The normal- 
ized variability measure is among the largest in the sample 
~ 17%/ 14% for the B/R bands (Table 1). As such, it is a 
promising candidate for photometric reverberation mapping. 

Figure \TT\ indicates that, at its redshift, the R band is rela- 
tively line-rich (the net contribution of the lines to the bands 
is ~ 6%), with Ha being the dominant emission line. We cal- 
culated <^ca an d <Haa using the ICCF and ZDCF algorithms. 
Unsurprisingly, the ZDC F method is so mewhat more noisy 
than the ICCF methods (Kas pi et alJl200(l see their Fig. 4). 
That said, we find the results of all methods of analyses to 
be consistent given the uncertainties. This means that inter- 
polation is not a major source of uncertainty for this object, 
which is consistent with the relative regular sampling of its 
light curves. 
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Figure 14. Photometric analysis of real B data and simulated R data for 
PG 0804. Several artificial time-lags (see legend) are assumed and the line- 
rich R band light curve is simulated using the line-poor B-band light curve 
(see §4.2.2). Analysis results for the real data are shown, for comparison, 
in blue shades. Both ICCF and ZDCF results are shown. Different artificial 
time-lags result in statistically different signals. As expected, longer time- 
lags result in a hump extending to longer timescales. The current data set is 
insensitive to short lags due to inadequate cadence. Qualitatively, an input 
time lag of 100-200 days is consistent with the observed signal. 

Our analysis shows the expected features: a hump in 
both ^ca an d <^aa> which peaks at similar (8t ~ 200 days) 
timescales. For comparison, we used a published light curve 
for the broad Ha emission line (iKaspi et alJ l^OOOFl and cal- 
culated %i c using the ZDCF algorithm. Clearly, both statisti- 
cal estimators defined here trace the line-to-continuum cross- 
correlation function with a varying degree of significance: 
^CA shows a > 2(7 peak (a being the uncertainty around the 
maximum), while the peak in ^aa is on ly ^ 1 o significant. 

Using the Monte Carlo simulations described in §4.2.1, 
peak or centroid time distributions were obtained, and are 
shown in figure [13] Clearly, time-lag estimates qualitatively 
agree between the different methods. Quantitatively, how- 
ever, there are some differences: the centroid method results 
in somewhat larger values for the time-lag. This results from 
the non-symmetric nature of E,ck: ^aa showing a more ex- 
tended wing toward longer time-delays (Fig. Q~2). Results 
based on the ZDCF approach result in errors, which are typ- 
ically smaller than those obtained using the ICCF method. 
This results from the fact that the ZDCF peak is largely de- 
termined by a few individual bins (note the two bins around 
200 days that lie above most other ZDCF points, as well as 
above the ICCF curve; Fig. [12). In fact, the resulting un- 
certainty on the ZDCF, being of order 20% in this case, is 
smaller than predicted based on our Monte Carlo simulations 
presented in figure [7] We conclude that estimating the un- 
certainty on the time-lag using the ZDCF results may be less 
reliable. 

Using the ICCF time-lag distributions, we determine the 
observed time-lag for PG0804 to be 250 ± 100 days, which 
is statistically consistent with the spectroscopically measured 
value of 193^ days (Kaspi et al. 2000). As far as systemat- 
ics is concerned, we do not expect a bias in the measured lag 
to exceed ~ 40% given that 'tot Admeasured £ 10 ( p ig- US- 
Using the numerical bootstrapping scheme defined in 
§4.2.2, we find that an input time lag of 100 — 200 days is qual- 
itatively consistent with the data. In particular, much shorter 
or much longer time-lags, are inconsistent with our findings in 

17 see http://wise-obs.tau.ac.il/~shai/PG/ 



terms of signal amplitude and shape. Interestingly, a weak yet 
marginally significant hump at around 200 — 300 days is im- 
plied by the data even when much shorter artificial time-lags 
are used as input. Therefore, one should be careful when in- 
terpreting marginally significant signals. For the specific case 
considered here we can conclude that the data cannot be used 
to detect short time-lags, which is of no surprise given that the 
cadence of the light curves is at 20 — 30 days intervals. 

To conclude, the time-lag estimated here by the photomet- 
ric reverberation mapping approach is ~ 200 ± 1 00 days, and 
is ther efore consistent with the value deduced by Kaspi et al. 
(2000), albeit with an expectedly larger uncertainty. 

4.3.2. PG 0844+349 

This low-z, low-luminosity object has a potentially promi- 
nent (~ 8%) line contribution to the /?-band from the Balmer 
emission lines. With < 70 points in its light curve, and what 
appears to be significant variations over 200 days time scales, 
our analysis reveals a marginal (> 1<7 using the FR-RSS 
method and > 2 using ZDCF statistics) peak at ~ 500 days 
(Fig. [151 ). The spectroscopically measured time-lag for this 
object is of order 50 days and agrees with the value expected 
from the BLR size vs. luminosity relation (table 1). Never- 
theless, the data do not show a corresponding peak at those 
time-scales. In fact, our simulations indicate that only a time- 
lag > 100 days, could have been detected with some cer- 
tainty. Interestingly, simulations with an input lag of 500 days 
(§4.2.2) are able to qualitatively account for the observed sig- 
nal (dashed magenta curve in figure [15). 

It is not clear whether a 500 days lag, ;/real, is associated 
with a long-term reverberation sign al of the BLR, sinc e there 
is no corresponding signal in the Kas pi et all (|2000) analy- 
sis of the spectroscopic data. We note, however, that it is 
also possible that the measurement uncertainties are some- 
what under-estimated for this object, hence the significance 
of our deduced signal over-estimated. That this might be the 
case is hinted by the normalized variability measure in both 
the B and R bands being similar (contrary to naive theoreti- 
cal expectations) and the fact that contamination by the host 
galaxy in this low-z quasar may be important (see below). 

4.3.3. PG 1229+204 

This quasar, being at low redshift (z = 0.064) has one of 
the most prominent (< 8%) Balmer line contributions to the 
/?-band among the PG objects in our sample. The analysis 
detects a significant trough at > 200 days, in both statistical 
estimators, and using both the ICCF and the ZDCF methods. 
This is unexpected since a peak rather than a trough is pre- 
dicted by the simulations for an expected lag of ~ 34 days 
(Fig. [15] ). Simulating cases with time-lags of order 200 days 
also produces a peak, in contrast to the observed signal. The 
interpretation that the B- rather than the /?-band is the line-rich 
band, is not supported by the o ptical spectrum sh owing rather 
typical quasar emission lines (Kaspi et al. 2000|). 

We cannot positively identify the origin of the unexpected 
signal but note that, at its low redshift, and being a rather faint 
low-luminosity quasar, the host galaxy is likely to substan- 
tially contribute to the flux in the bands. This has a dual result: 
1) the actual contribution of the emission lines to the bands is 
smaller than that predicted here, and 2) the data may be sub- 
stantially affected by varying seeing conditions between the 
nights, thereby contributing to the variance in the light curves, 
and masking the true signal. Interestingly, the reduced vari- 
ability measure in both bands is similar for this object, which 
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Figure 15. Analysis of 15 PG quasars from the sample of Kaspi et al. 1 2000, see our table 1). In each pair of columns, the left panel shows the B (blue points) 
and R (red points) light curves, and the right panel the calculated £ca (blue curves and shades) and ^aa (red curves and shades) using the ICCF (solid lines) 
and ZDCF (points with error-bars) numerical schemes. Cases for which a statistically significant signal is detected in both estimators and using both numerical 
schemes, are highlighted. For each quasar, the results of three £ca simulations are shown (in dashed blue-shaded curves) based on the expected time-lag from the 
BLR size vs. luminosity relation, (j' ag ex (table 1): light to dark shades correspond simulated time-lags of cx /2 — > fj^ cx — > 2f' ag cx (cyan curves correspond to 
'lag cx)- F° r PG0844, a simulation with a 500 days lag is also included (dashed magenta line). The ICCF results reported here are used to construct the ensemble 
averages considered in figure H71 
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may indicate the presence of an additional source of noise not 
accounted for by the reported measurement uncertainties. 

4.3.4. PG 1613+658 

Analyzing the data for this object, a marginally significant 
signal (at the 1 a — 2a level) is detected in both statistical es- 
timators, at < 300 days (Fig. [TBI ). Formally, £,qa statistics 
implies that the time lag is 300 ± 130 days. At its redshift, 
the net contribution of emission lines is of order 5%, and is 
therefore larger than the typical flux measurement uncertainty 
in this object (< 2%). Simulations with an artificial lag of 
~ 300 days are able to qualitatively account for the observed 
signal. Much shorter time lags [e . g., ~ 40 days, as spectro- 
scopically deduced bv lKaspi et alj ( 2000)] do not seem to be 
supported by our findings. Interestingly, for this object, spec- 
troscopic time-lag measurements using the peak distribution 
function give very large errors for the Ha line, and given its 
luminosity, PG 1613+658 l ies considerably below the BLR 
size vs. luminosity relation dKaspi et al.ll2.00Ql see also our ta- 
ble 1). Using our deduced lag, a better agreement is obtained 
with the Kaspi et all (120001) size-luminosity relation. 

4.3.5. PG 1700+518 

This object shows a statistically significant (2a — 3a), wide 
peak extending over the range 150-450days. Assuming that 
Hj3 is the main line contributor to the R band's flux, our simu- 
lations indicate that the signal is consistent with a time-lag of 
~ 430 days (formally, £,ca statistics yields a 300 ± 160days' 
lag). This value is statistically cons istent with the time-delay 
measurement of iKaspi et alj £2000), but our deduced lag is 
in somewhat better agreement with the expected BLR size 
given the source luminosity (table 1). Our simulations show 
that time-lags considerably shorter (<C 200 days) or longer 
(S> 800 days) are less favored by the data and cannot fully 
account for the observed signal (Fig. H3l . 

4.3.6. PG 1704+608 

This object, having ~ 40 visits in each band, shows a 
marginal (> la) hump at 200-300 days. Formally, the de- 
duced time -lag is 400 ± 200 days, and is consistent with the 
values of Kas pi et all d2000h . At z = 0.371, the net contri- 
bution of emission lines to the line-rich 7?-band is only about 
twice as large as the mean measurement uncertainty. This, 
and the small number of data points, account for the low- 
significance of the observed signal. Further, simulations show 
that a marginal signal is indeed expected to occur at roughly 
the spectroscopically determined lag. 

4.3.7. PG 21 30+099 

There is a hint for a localized (~ la) peak around 200- 
300 days in both statistical estimators (£,ca statistics yields a 
lag of 240 ± 100 days) using both the ICCF and the ZDCF 
methods. At face value, this timescale is consistent with 
the results of IKaspi et al.l (120001) who fi nd a time lag of ~ 
200 days. Nevert heless, as discusse d by iGrier et al.N 2008). 
the cadence of the Kaspi et al. (2000) observations may not be 
adequ ate to uncover the true time lag in this object. Specifi- 
cally, |GrieretaD (EQQD) find a time lag of order 30 days using 
better sampled light curves on short time scales. The fact that 
we do not find a significant peak on short ti mescales is not sur- 
prising since we are essentially using the IKaspi et ail d2000l) 
data set. 



4.3.8. PG0026+129 and other objects yielding null results 

Like PG 0804, PG 0026+129 (z = 0. 142) is among the best- 
sampled quasars in our sample with ~ 70 points per light 
curve (see Fig. [TBI . While the normalized variability mea- 
sure is comparable to that of PG 0804, most of the power is 
associated with long-term variability. In particular, when sub- 
tracting a linear trend from the light curve, the normalized 
variability measure drops by a factor of ~ 2 in both bands. At 
its redshift, the /?-band is the line-rich (Fig. QT} with a net 
relative line contribution of ~ 3% due to Ha (or ~ 7% if iron 
is important). 

The statistical estimators calculated here do not yield a 
significant time-lag. Specifically, there is no discernible 
peak detected in either estimator using either cross-correlation 
method (i.e., ICCF or ZDCF). While de-trending seems to 
somewhat improve the signal (there is a 1 a hump at around 
100 days), it falls short of revealing a significant result. We 
attribute our failure to securely identify a signal to a combi- 
nation of factors related to the smaller variability measure on 
the relevant BLR-size timescales, as well as to the lower net 
contribution of emission lines to the filters, compared to the 
case of PG 0804. 

PG 0052+251 has 72 points in each band, and a normal- 
ized variability measure of 26% /20% in the B/R bands, partly 
due to a sharp feature extending over ~ 200 days scales. The 
typical measurement uncertainty is < 3% in both bands. At 
z = 0.155, the net line contribution is of order 4%. i.e., com- 
parable to the noise-level. Therefore, despite the relatively 
large number of visits, no clear signal is detected, in either 
statistical estimator, using either numerical s c heme . 

The remaining quasars in the IKaspi et all d2000h PG sam- 
ple also do not yield a significant signal. This is however 
expected given the small number of points in the photomet- 
ric light curve of these quasars being ~40 (see Fig. |7). The 
fact that such objects do have a spectroscopically-determined 
time lag is also not surprising: the photometric approach to 
reverberation mapping is less sensitive than the spectroscopic 
one. 

4.3.9. Ensemble averages 

To overcome the limitations plaguing the analysis of in- 
dividual objects, ensemble averages may be considered. As 
demonstrated in §3, combining the results for many objects, 
allows the time-lag to be determined with good accuracy, and 
was a main motivation on our part to develop this approach to 
reverberation mapping. Clearly, large ensembles are required 
to reach a robust result when less than favorable light curves 
are available on an individual case basis. Nevertheless, with 
upcoming photometric surveys, statistics is unlikely to be a 
limiting factor, and defining ensemble averages over narrow 
luminosity and redshift bins will be possible. 

As our small sample of 17 PG quasars consists of a range 
of redshifts and quasar luminosities, we use the follow- 
ing transformation to put the results for different objects 
on an equal footing: 8t — > St /[(I +z)Lj 5 ], where L45 = 

XLx (5 1 00 A) / 1 45 erg s ~ 1 . This transformation corrects for 
both cosmological time-dilation, an d the fact th a t the BLR 
size scales roughly as Z/. Following IKaspi et all (120001) . we 
take 7 = 0.7 but note that sim ilar results are obtained for 
0-5 < 7 < 0.7 dBentz et alj|2009l and references therein). Us- 
ing the above transformation, all PG quasars are effectively 
transformed to a z — quasar with L45 = 1 . 
Figure [17] shows the ensemble-averaged results for 
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Figure 16. Light curve analysis for PG 0026+129. Left: The B and R light curves show significant variability on timescales of the order the duration of the light 
curve. De-trended light curves were calculated by subtracting a l sl -degree polynomial. The normalized variability measure for the de-trended light curves is a 
factor two smaller than the orig inal data. Right: the calculated statistical estimators for the de-trended light curves (upper panel) and the original data (lower 
panel). See the caption of figure [T2l for the meaning of different curves. While de-trending seems to be able to improve the signal, a statistically significant peak 
is not detected by either method. 

occurring at times which are qualitatively co nsistent with the 
BLR sizes implied by the Kas pi et ail ([2000) results. A larger 
sample of objects will allow to better quantify the mean time 
lag for sub-samples of quasars over narrow luminosity and 
redshift intervals. Alternatively, better quality data (higher 
S/N, and better sampled light curves) could be used to better 
constrain the time-lag for individual objects in the PG sample 
of quasars. 

5. SUMMARY 

We demonstrated that line-to-continuum time-delays asso- 
ciated with the broad emission line region in quasars can be 
deduced from broadband photometric light curves. This is 
made possible by defining appropriate statistical estimators 
[termed 4ca(5?) and ^AA(5f)] to process the data, and which 
effectively subtract the continuum contribution to the cross- 
correlation signal of the light curves. Using this formalism, 
the time-lag is identified with the time 8t at which these esti- 
mators peak. 

The proposed formalism makes use of the fact that line and 
continuum processes have very different variability proper- 
ties, and that quasars exhibit large amplitude variations on 
timescales comparable to the BLR light crossing time. The 
formalism is generally adequate for cases in which few strong 
emission lines, having different intensities, contribute to the 
spectrum, as indeed characterizes the optical emission spec- 
trum of quasars. 

To effectively apply the method, prior knowledge of the 
quasar redshift (using either photometric or spectroscopic 
means) is required to identify relatively line-dominated and 
line-free spectral bands. Together with knowledge of the fil- 
ter response functions, these may be used to better implement 
the algorithm and interpret the results. Using numerical simu- 
lations, we demonstrated that photometric reverberation map- 
ping is feasible under realistic observing conditions and for 
typical quasar properties. We showed, by means of numer- 
ical simulations and real data (pertaining to the PG quasar 
sample), that the measurement of the line-to-continuum time- 
delay is robust although somewhat less accurate than that ob- 
tained using the spectroscopic reverberation mapping tech- 
nique. In addition, we quantify biases in the photometrically 
deduced time-lags, which are sensitive to the BLR geometry 
and the duration of the light curve. We identify their origin 
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Figure 17. Ensemble averages for £,ch (upper panel) and i^ aa (lower panel), 
for the PG sample of 17 quasars from Kaspi et al. 12000, see our table 1). 
Both the arithmetic mean (dark shades) and the median (light shades) are 
shown. Ensemble averages for both statistical estimators give consistent re- 
sults, and peak around 180-280 days. This is in agreeme nt with the predicte d 
time-lag based on the BLR size vs. luminosity relation of Kaspi et al. 1 2000), 
for a typical L45 = 1 quasar (see text). The median shows a less significant 
peak although the results are qualitatively consistent with the ensemble mean. 
This demonstrates the feasibility of the photometric reverberation mapping 
technique in uncovering the typical time-lag in an ensemble of objects for 
which individual time-lag determination is difficult to achieve. 

fcAi <^aa f° r the entire PG sample (table 1). We deliber- 
ately include all objects regardless of data quality issues (as 
in the case of e.g., PG 1229). A peak in both (mean) statisti- 
cal estimators is seen around 200 days. This value is similar 
to a time-lag of 17 days, which is the value predicted by the 
iKasoi et alJ (12000) BLR size-luminosity relation for a L45 = 1 
quasar. Also plotted are the median results for both statistical 
estimators. Both the mean and the median provide qualita- 
tively consistent results, although quantitative differences do 
exist. Given the small nature of the sample, we do not present 
estimates for the uncertainty on the mean/median values. 

To conclude, averaging the results for the PG quasar sam- 
ple, we obtain a similar peak in both statistical estimators, 
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and suggest several methods to correct for them. In addi- 
tion, we qualitatively study potential biases associated with 
the presence of multiple emission lines in quasar spectra, and 
with finite inter-band continuum time-delays. 

The main advantage of the photometric approach to rever- 
beration mapping over other methods is that it is consider- 
ably more efficient and cheap to implement, and can there- 
fore be applied to large samples of objects. In fact, any sur- 
vey which photometrically monitors the sky with proper ca- 
dence and signal-to-noise in two or more filters, and for which 
quasars may be identified and their redshift estimated, can be 
used to measure the size of the BLR. We demonstrated that, 
with large enough samples of quasars, it is possible to beat 
down the noise, and deduce (statistical) BLR size measure- 
ments that are considerably more accurate than are currently 
achievable by spectroscopic means. 

By combining photometric light curves with single epoch 
spectroscopic measurements, which allow to estimate the ve- 
locity dispersion of the BLR, it may be possible to measure 
the black hole mass in orders of magnitude more objects than 
are directly measurable by other means. This would allow 
for the statistical measurement of the SMBH mass for sub- 
samples of quasars (selected according to e.g., their luminos- 
ity, redshift, or colors), leading to highly accurate results, and 
potentially revealing differences between different classes of 
objects, which, thus far, may have eluded detection. 

To conclude, photometric surveys that monitor the sky on 
a regular basis are likely to play a major role in the coming 
decades in the (statistical) measurement of the BLR size in 
quasars, as well as in determining their black hole masses. 
This will lead to a more complete census of black holes over 
cosmic time, and shed light on outstanding problems concern- 
ing the formation and co-evolution of black holes and their 
host galaxies. 
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